home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_09_05 / 9n05085a < prev    next >
Text File  |  1991-03-20  |  12KB  |  269 lines
  1.  
  2. /**********************************************************
  3.  *
  4.  * Copyright (c) 1991, John Forkosh.  All rights reserved.
  5.  * Released to the public domain only for use that is both
  6.  * (a) by an individual, and (b) not for profit.
  7.  * --------------------------------------------------------
  8.  *
  9.  * Function:    lint1d ( x0,dx,n, f,y )
  10.  *
  11.  * Purpose:     Creates a table containing n values,
  12.  *              starting with y[0]=f(x0), y[1]=f(x0+dx),
  13.  *              through y[n-1]=f(x0+(n-1)*dx).
  14.  *              However, the value in y[i] won't
  15.  *              necessarily exactly equal f(x0+i*dx)
  16.  *              as indicated above.  Instead, the table
  17.  *              is constucted to minimize the mean square
  18.  *              error due to sampling random points in
  19.  *              f(x)'s domain.
  20.  *
  21.  * Arguments:
  22.  *              x0 (I)  Double containing the first point
  23.  *                      in the function's domain.
  24.  *              dx (I)  Double containing x-interval
  25.  *                      between points in the table.
  26.  *              n (I)   Int containing the number of
  27.  *                      points in the table.
  28.  *              f (I)   Address of function returning
  29.  *                      double whose values are to be
  30.  *                      represented in the table.
  31.  *              y (O)   Address of double returning the
  32.  *                      table of n values, as discussed.
  33.  *
  34.  * Returns:     (int)   0 for a normal solution, or
  35.  *                      1 for any error (e.g., a singular
  36.  *                      set of equations).
  37.  *
  38.  * Notes:     o See the accompanying article for a complete
  39.  *              description of lint1d().
  40.  *            o Set #defined value of TESTDRIVE to 1/TRUE
  41.  *              to compile a sample main() (see below).
  42.  *
  43.  * Source:      LINT1D.C
  44.  *
  45.  * --------------------------------------------------------
  46.  * Revision History:
  47.  *
  48.  * 01/07/91     J.Forkosh       Installation
  49.  *
  50.  *********************************************************/
  51.  
  52. /* --- standard headers --- */
  53. #include <stdio.h>              /* need NULL ptr value */
  54. #include <stdlib.h>             /*for malloc() prototype*/
  55. #define msglevel 1              /* to printf debug info */
  56.  
  57. /* --- set testdrive to compile (or not) test main() --- */
  58. #define TESTDRIVE 1             /* 1=compile it (0=not) */
  59.  
  60. /* --------------------------------------------------------
  61. Entry point
  62. -------------------------------------------------------- */
  63. int     lint1d ( x0,dx,n, f,y ) /* returns 0=ok, 1=error */
  64. double  x0,dx;                  /*1st point and interval*/
  65. int     n;                      /* number of points */
  66. double  (*f)();                 /*func to be interpolated*/
  67. double  *y;                     /*table for interpolation*/
  68. {
  69. /* --- local allocations and declarations --- */
  70. int     iserror = 0;            /* no error detected yet */
  71. int     i,j;                    /* row,col indexes */
  72. double  x;                      /* arg for f(x) */
  73. /* --- need temporary array for coefficient matrix --- */
  74. double  *a = (double *)malloc((n*n+1)*sizeof(double));
  75.  
  76. /* --------------------------------------------------------
  77. Set up coefficient matrix as per discussion in article
  78. (since the matrix is symmetric, the columnwise requirement
  79. for simq() input is irrelevant).
  80. -------------------------------------------------------- */
  81. /* --- first check that malloc() was successful --- */
  82. if ( a == NULL )                /*failed to malloc matrix*/
  83.   return ( iserror=1 );         /*return error to caller*/
  84. for (i=0;i<n*n;i++) a[i]=0.0;   /* init array to zeros */
  85. /* --- 4's on diag (2's at extremes) and 1's offdiag --- */
  86. for ( i=1; i<n; i++ )           /* skip 0,0 element */
  87.   {
  88.   j = n*i + i;                  /* index of diag element */
  89.   a[j] = 4.0;                   /* set diagonal element */
  90.   a[j-1] = a[j+1] = 1.0;        /* and off-diag elements */
  91.   } /* --- end-of-for(i=1...n-1) --- */
  92. a[0] = a[n*n-1] = 2.0;          /* set 1st,last diag */
  93. a[1] = a[n*n-2] = 1.0;          /*and remaining off-diags*/
  94.  
  95. /* --------------------------------------------------------
  96. Set up vector of constants as per discussion in article
  97. -------------------------------------------------------- */
  98. /* --- interior points --- */
  99. for ( i=1; i<n; i++ )           /* loop skips 1st point */
  100.   {
  101.   x = x0 + dx*i;                /* next arg to be tabled */
  102.   y[i] = 2.0*(f(x-0.5*dx)+      /*const vector calculated*/
  103.                 f(x)+f(x+0.5*dx)); /* as per article */
  104.   } /* --- end-of-for(i=1...n-1) --- */
  105. /* --- boundary points --- */
  106. y[0]   = f(x0) + 2.0*f(x0+0.5*dx); /* first x in domain */
  107. y[n-1] = f(x)  + 2.0*f(x -0.5*dx); /*use last x from loop*/
  108.  
  109. /* --------------------------------------------------------
  110. Display input to simq (for debugging purposes if necessary)
  111. -------------------------------------------------------- */
  112. if ( msglevel >= 9              /* want debugging output */
  113. &&   n < 10 )                   /* have room to fit it */
  114.   {
  115.   printf("lint1d> input to simq...\n"); /* stub info */
  116.   for ( i=0; i<n; i++ )
  117.     { /* --- display row (really col) of the matrix --- */
  118.       for ( j=0; j<n; j++ ) printf("%6.2lf",a[n*i+j]);
  119.       /* --- followed by literal for y* and const y --- */
  120.       printf("    y*[%d]    %8.2lf\n", i,y[i]); }
  121.   } /* --- end-of-if(msglevel>=9) --- */
  122.  
  123. /* --------------------------------------------------------
  124. Solve the simultaneous equations
  125. -------------------------------------------------------- */
  126. iserror = simq(a,y,n);          /* y[] returns solution */
  127. free ( (void *)a );             /*don't need coeff matrix*/
  128. return ( iserror );             /*back to caller with y[]*/
  129. } /* --- end-of-function lint1d() --- */
  130.  
  131.  
  132. #if TESTDRIVE
  133. /**********************************************************
  134.  *
  135.  * Purpose:     Test driver for lint1d().
  136.  *              Prepares a table of values for
  137.  *              interpolation of the function f(x)=x*x,
  138.  *              and compares the resulting accuracy
  139.  *              with usual linear interpolation.
  140.  *
  141.  * Command-Line Arguments:
  142.  *              x0 (I)  Double containing the first point
  143.  *                      in the function's domain
  144.  *                      (default=-10.0).
  145.  *              dx (I)  Double containing x-interval
  146.  *                      between points in the table
  147.  *                      (default=1.0).
  148.  *              n (I)   Int containing the number of
  149.  *                      points in the table
  150.  *                      (default=21).
  151.  *              expon (I) Int containing the exponent for
  152.  *                      the test function f(x)=x**expon
  153.  *                      (default=2).
  154.  *
  155.  *********************************************************/
  156.  
  157. /* --- program defaults (x0,dx,n,expon from cmdline) --- */
  158. #define X0 (-10.0)              /* 1st point in table */
  159. #define DX 1.0                  /* table interval */
  160. #define N 21                    /*number pts in table*/
  161. static  int expon=2;            /* test f(x)=x**expon */
  162. #define NMAX 99                 /* easier than malloc */
  163. #define NRMS 101                /*pts per dx for rms calc*/
  164.  
  165. /* --------------------------------------------------------
  166. Entry point
  167. -------------------------------------------------------- */
  168. main    ( argc, argv )
  169. int     argc;
  170. char    *argv[];
  171. {
  172. /* --- local allocations and declarations --- */
  173. double  x0=X0, dx=DX; int n=N;  /* defaults */
  174. int     i,j;                    /* loop indexes */
  175. double  x;                      /* arg for f() */
  176. double  f();                    /* test function */
  177. double  y[NMAX];                /* table from lint1d() */
  178. int     iserror=0;              /*return flag from lint1d*/
  179. double  frms=0.0,yrms=0.0;      /* mean square errors */
  180. double  xa,xb;                  /*interval bounds for rms*/
  181. double  sqrarg;                 /*dummy arg for sqr macro*/
  182. /* --- linear interpolation (u(x=xa)=ya, u(x=xb)=yb) --- */
  183. #define u(x,xa,ya,xb,yb) (ya*(xb-x)+yb*(x-xa))/(xb-xa)
  184. /* --- square a double --- */
  185. #define sqr(x) (sqrarg=(x),sqrarg*sqrarg)
  186.  
  187. /